General project set-up

# Load libraries and sources required to run the script
    library(tidyverse)
    library(ggthemes)
    library(lmerTest)
    library(emmeans)
    library(multcomp)
    library(effects)
    library(gridExtra)
    library(rstatix)

# Default ggplot settings

    Fill.colour<-scale_colour_manual(values = c("black", "gray70", "gray35"))

    ggthe_bw<-theme(plot.background=element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          legend.box.background = element_rect(),
          panel.background =element_rect(fill = NA, color = "black")
          )+
    theme_bw()

Data exploration

1. Get the files with all the YII by species

  YII.data<-read.csv("YII_Data/All_YII_data.csv", header = T)
  summary(YII.data)
##         Sample             Date      Spp          Fragment    Treatment 
##  Ac_288_T21:   2   2017-11-16: 354   Ac:2281   Ac_108 :  24   A  :2096  
##  Ac_101_T0 :   1   2017-11-23: 354   Of:1051   Ac_116 :  24   N  :1974  
##  Ac_101_T1 :   1   2017-11-29: 354   Ss:2687   Ac_119 :  24   N+P:1949  
##  Ac_101_T10:   1   2017-12-06: 354             Ac_122 :  24             
##  Ac_101_T11:   1   2017-12-13: 354             Ac_143 :  24             
##  Ac_101_T12:   1   2018-01-03: 354             Ac_152 :  24             
##  (Other)   :6012   (Other)   :3895             (Other):5875             
##  Replicate      YII            Genotype         Days        
##  R1:3061   Min.   :0.0000   G_62   : 577   Min.   :-112.00  
##  R2:2958   1st Qu.:0.4420   G_48   : 570   1st Qu.:  14.00  
##            Median :0.5030   G_07   : 462   Median :  65.00  
##            Mean   :0.4988   Ss_20  : 402   Mean   :  52.42  
##            3rd Qu.:0.5720   Ss_23  : 402   3rd Qu.:  92.00  
##            Max.   :0.6870   Ss_27  : 400   Max.   : 156.00  
##                             (Other):3206                    
##    Time_Point         Phase         TotalSH          logSH       
##  T10    : 354   Baseline : 954   Min.   :0.000   Min.   :-5.473  
##  T5     : 354   Heat     :1520   1st Qu.:0.020   1st Qu.:-1.709  
##  T6     : 354   Nutrients:2818   Median :0.057   Median :-1.243  
##  T7     : 354   Ramping  : 514   Mean   :0.116   Mean   :-1.317  
##  T8     : 354   Recovery : 213   3rd Qu.:0.158   3rd Qu.:-0.802  
##  T9     : 354                    Max.   :1.394   Max.   : 0.144  
##  (Other):3895                    NA's   :4657    NA's   :4657    
##      D.Prp        Community InitialCommunity
##  Min.   :0.0000   A :2281   A :2281         
##  1st Qu.:0.0000   B : 195   B : 186         
##  Median :0.0000   C1: 538   C1: 468         
##  Mean   :0.1705   C3: 727   C3: 762         
##  3rd Qu.:0.0000   D :2278   D :2322         
##  Max.   :1.0000                             
##  NA's   :2836

Merge/Transform

# Organize data type
      YII.data$Date<-as.Date(YII.data$Date, "%Y-%m-%d")
      YII.data$Days<-(as.numeric(YII.data$Date) -17485)
      #Time as a factor, not as int
      str(YII.data)
## 'data.frame':    6019 obs. of  16 variables:
##  $ Sample          : Factor w/ 6018 levels "Ac_101_T0","Ac_101_T1",..: 1 2 10 11 12 13 14 15 16 17 ...
##  $ Date            : Date, format: "2017-07-26" "2017-08-30" ...
##  $ Spp             : Factor w/ 3 levels "Ac","Of","Ss": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Fragment        : Factor w/ 354 levels "Ac_101","Ac_102",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment       : Factor w/ 3 levels "A","N","N+P": 3 3 3 3 3 3 3 3 3 3 ...
##  $ Replicate       : Factor w/ 2 levels "R1","R2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ YII             : num  0.644 0.576 0.563 0.568 0.645 0.589 0.595 0.606 0.605 0.606 ...
##  $ Genotype        : Factor w/ 17 levels "G_07","G_08",..: 5 5 5 5 5 5 5 5 5 5 ...
##  $ Days            : num  -112 -77 -36 -27 -9 1 8 14 21 28 ...
##  $ Time_Point      : Factor w/ 25 levels "T0","T1","T10",..: 1 2 12 19 20 21 22 23 24 25 ...
##  $ Phase           : Factor w/ 5 levels "Baseline","Heat",..: 1 1 1 1 1 1 3 3 3 3 ...
##  $ TotalSH         : num  NA NA NA NA NA ...
##  $ logSH           : num  NA NA NA NA NA ...
##  $ D.Prp           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Community       : Factor w/ 5 levels "A","B","C1","C3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ InitialCommunity: Factor w/ 5 levels "A","B","C1","C3",..: 1 1 1 1 1 1 1 1 1 1 ...
      YII.data$DaysF<-as.factor(YII.data$Days)

      
      YII.data$Spp <- as.factor(YII.data$Spp)
      
      YII.data$Treatment <- as.factor(YII.data$Treatment)
      
      YII.data$Genotype<-factor(as.character(YII.data$Genotype), 
                             levels=c("G_48", "G_62","G_31", 
                                      "G_08","G_07", "G_50", 
                            "Of_34","Of_20","Of_6", "Of_31",
                            "Ss_22","Ss_23","Ss_27", "Ss_28",
                            "Ss_20", "Ss_24","Ss_30"
                             ))  # D dominance order
      YII.data$Community <- factor(YII.data$Community, 
                                   levels = c("A", "B", "C3", "C1", "D"))
      
# Check the data
      str(YII.data)
## 'data.frame':    6019 obs. of  17 variables:
##  $ Sample          : Factor w/ 6018 levels "Ac_101_T0","Ac_101_T1",..: 1 2 10 11 12 13 14 15 16 17 ...
##  $ Date            : Date, format: "2017-07-26" "2017-08-30" ...
##  $ Spp             : Factor w/ 3 levels "Ac","Of","Ss": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Fragment        : Factor w/ 354 levels "Ac_101","Ac_102",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment       : Factor w/ 3 levels "A","N","N+P": 3 3 3 3 3 3 3 3 3 3 ...
##  $ Replicate       : Factor w/ 2 levels "R1","R2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ YII             : num  0.644 0.576 0.563 0.568 0.645 0.589 0.595 0.606 0.605 0.606 ...
##  $ Genotype        : Factor w/ 17 levels "G_48","G_62",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ Days            : num  -112 -77 -36 -27 -9 1 8 14 21 28 ...
##  $ Time_Point      : Factor w/ 25 levels "T0","T1","T10",..: 1 2 12 19 20 21 22 23 24 25 ...
##  $ Phase           : Factor w/ 5 levels "Baseline","Heat",..: 1 1 1 1 1 1 3 3 3 3 ...
##  $ TotalSH         : num  NA NA NA NA NA ...
##  $ logSH           : num  NA NA NA NA NA ...
##  $ D.Prp           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Community       : Factor w/ 5 levels "A","B","C3","C1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ InitialCommunity: Factor w/ 5 levels "A","B","C1","C3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ DaysF           : Factor w/ 26 levels "-112","-77","-36",..: 1 2 3 4 5 6 7 8 9 10 ...
      summary(YII.data)
##         Sample          Date            Spp          Fragment   
##  Ac_288_T21:   2   Min.   :2017-07-26   Ac:2281   Ac_108 :  24  
##  Ac_101_T0 :   1   1st Qu.:2017-11-29   Of:1051   Ac_116 :  24  
##  Ac_101_T1 :   1   Median :2018-01-19   Ss:2687   Ac_119 :  24  
##  Ac_101_T10:   1   Mean   :2018-01-06             Ac_122 :  24  
##  Ac_101_T11:   1   3rd Qu.:2018-02-15             Ac_143 :  24  
##  Ac_101_T12:   1   Max.   :2018-04-20             Ac_152 :  24  
##  (Other)   :6012                                  (Other):5875  
##  Treatment  Replicate      YII            Genotype         Days        
##  A  :2096   R1:3061   Min.   :0.0000   G_62   : 577   Min.   :-112.00  
##  N  :1974   R2:2958   1st Qu.:0.4420   G_48   : 570   1st Qu.:  14.00  
##  N+P:1949             Median :0.5030   G_07   : 462   Median :  65.00  
##                       Mean   :0.4988   Ss_23  : 402   Mean   :  52.42  
##                       3rd Qu.:0.5720   Ss_20  : 402   3rd Qu.:  92.00  
##                       Max.   :0.6870   Ss_27  : 400   Max.   : 156.00  
##                                        (Other):3206                    
##    Time_Point         Phase         TotalSH          logSH       
##  T10    : 354   Baseline : 954   Min.   :0.000   Min.   :-5.473  
##  T5     : 354   Heat     :1520   1st Qu.:0.020   1st Qu.:-1.709  
##  T6     : 354   Nutrients:2818   Median :0.057   Median :-1.243  
##  T7     : 354   Ramping  : 514   Mean   :0.116   Mean   :-1.317  
##  T8     : 354   Recovery : 213   3rd Qu.:0.158   3rd Qu.:-0.802  
##  T9     : 354                    Max.   :1.394   Max.   : 0.144  
##  (Other):3895                    NA's   :4657    NA's   :4657    
##      D.Prp        Community InitialCommunity     DaysF     
##  Min.   :0.0000   A :2281   A :2281          1      : 354  
##  1st Qu.:0.0000   B : 195   B : 186          8      : 354  
##  Median :0.0000   C3: 727   C1: 468          14     : 354  
##  Mean   :0.1705   C1: 538   C3: 762          21     : 354  
##  3rd Qu.:0.0000   D :2278   D :2322          28     : 354  
##  Max.   :1.0000                              49     : 354  
##  NA's   :2836                                (Other):3895

Remove / subset timepoints

 # Remove baseline values
      YII.data<-subset(YII.data, Days>-1)
    # Remove recovery values
      YII.data<-subset(YII.data, Days<112)
      # write.csv(YII.data, "Outputs/Experiment_YII_data.csv", row.names = F)  
    # YII.Wide<- reshape(YII.data, idvar = "Fragment", timevar = "Days", direction = "wide")
      
      Spp.fragments<-YII.data %>% 
        group_by(Spp, Genotype, Treatment, Replicate) %>% count(Fragment)
      Spp.fragments
      #write.csv(Spp.fragments, "Outputs/Meassurments_perFragments.csv", row.names = F)

# Subset data 
      YII.nutrients<-subset(YII.data, Days<80)
      YII.heat<-subset(YII.data, Days>75)

2. Exploratory graphs

All time points (nutrients + heat stress)

  • Colony (Genotype) differences
YII_Colony<- ggplot(YII.data, aes (Days, YII, colour=Genotype)) +
        stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.5)+
        stat_summary(fun.y=mean, geom="line", alpha=0.6) +   theme_bw()
      YII_Colony + ylim(0.0, 0.65) + facet_grid (Spp~Treatment)

YII_Frag_Gen<- ggplot(YII.data, aes (Days, YII, colour=Genotype, group=Fragment)) +
        stat_summary(fun.y=mean, geom="line", alpha=0.5) +  
        theme_bw() + theme(legend.position = "bottom",
                           legend.title = element_blank())
      YII_Frag_Gen + ylim(0.0, 0.65) + facet_grid (Spp~Treatment)

Figure 3: Treatment differences by species and dominant symbiont

DropPlot<-YII.data[which(YII.data$Spp !="Ac" &
                        YII.data$DaysF==110),]

YII.datab<-YII.data[!(YII.data$Sample) %in% (DropPlot$Sample),]
  
YII_Treat<- ggplot(YII.datab, aes (Days, YII, colour=Treatment)) +
        #stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
        #             width = 0.2, position = position_dodge(1) ) +
        ggthe_bw+ Fill.colour +
        # stat_summary(fun.y=mean, geom="point")  +
        theme(legend.position="bottom",
              strip.background = element_rect(fill="white"))+
        scale_y_continuous(limits = c(0, 0.7),
                           breaks = seq(0, 0.6, 0.1),  
                           expand = c(0, 0),
                           name=("Fv/Fm")) +
        scale_x_continuous(name="Days in the experiment",
                           limits = c(-1,113),
                           breaks = seq(0, 113, 15),  
                           expand = c(0, 0))+
        annotate("segment", x = 2, xend = 91, y = 0.05, yend = 0.05,
                  colour = "gray90", linetype=2)+
        annotate("segment", x = 79, xend = 91, y = 0.05, yend = 0.20,
                  colour = "gray90", linetype=3)+
        annotate("segment", x = 91, xend = 110, y = 0.20, yend = 0.20,
                  colour = "gray90", linetype=3)+
        annotate("text", x = 45, y = 0.12, label = "Nutrients", size=3)+
        annotate("text", x = 99, y = 0.12, label = "Heat", size=3)
#YII_Treat + facet_grid (~Spp) + geom_smooth(span=0.5)

All_SppFif<-YII_Treat + facet_wrap (Spp~Community) + geom_smooth(span=0.5)
All_SppFif

#ggsave(file="Outputs/All_spp_YII_Treatb.svg", plot=All_SppFif, width=7, height=6)

Figure 3b (No smooth)

YII_Treat_BW<- ggplot(data=YII.data, aes (Days, YII, colour=factor(Treatment), shape=factor(Treatment))) + 
        ggthe_bw + Fill.colour+
        stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                     position = position_dodge(1) )+
        stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                     linetype=1, alpha=1) + 
         stat_summary(fun.y=mean, geom="point", size = 2,
                   position=position_dodge(width=1), alpha=0.8)  +
        theme(legend.position=c(0.1, 0.8),
        legend.title = element_blank(), 
        strip.background =element_rect(fill=NA)) + # geom_smooth()+
      
    scale_y_continuous(limits = c(0.1, 0.7),
                           breaks = seq(0.1, 0.6, 0.1),  
                           expand = c(0, 0),
                           name=expression(~italic("Fv / Fm"))) +
        scale_x_continuous(name="Days in the experiment",
                           limits = c(-1,113),
                           breaks = seq(0, 113, 15),  
                           expand = c(0, 0))+
        annotate("segment", x = 2, xend = 91, y = 0.12, yend = 0.12,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 79, xend = 91, y = 0.12, yend = 0.20,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 91, xend = 110, y = 0.20, yend = 0.20,
                  colour = "gray90", linetype=1)
        
      
  Figure3b<-YII_Treat_BW + #facet_grid (Spp~.)+
        annotate("text", x = 45, y = 0.15, label = "Nutrients", size=3)+
        annotate("text", x = 99, y = 0.15, label = "Heat", size=3)+
        facet_wrap(Spp~Community)
  Figure3b

#ggsave(file="Outputs/Fig_3b_Acer_YII_Treat.svg", plot=Figure3, width=3.5, height=3)

Figure 3c (Symbiont based)

YII_Treat_Sym_Spp<- ggplot(data=subset(YII.data), aes (Days, YII, colour=factor(InitialCommunity))) +
         ggthe_bw + geom_smooth()+
   theme(legend.position="bottom", 
              legend.title = element_blank(),
              strip.background = element_rect(fill="white"))+
  scale_y_continuous(limits = c(0.2, 0.7),
                           breaks = seq(0.1, 0.65, 0.1),  
                           expand = c(0, 0),
                           name=expression(~italic("Fv / Fm"))) +
  scale_x_continuous(name="Days in the experiment",
                           limits = c(-1,113),
                           breaks = seq(0, 113, 15),  
                           expand = c(0, 0))+
  annotate("segment", x = 2, xend = 91, y = 0.21, yend = 0.21,
                  colour = "gray90", linetype=2)+
  annotate("segment", x = 79, xend = 91, y = 0.21, yend = 0.28,
                  colour = "gray90", linetype=3)+
  annotate("segment", x = 91, xend = 110, y = 0.28, yend = 0.28,
                  colour = "gray90", linetype=3)

        
Figure3c<-YII_Treat_Sym_Spp + facet_grid (Spp~Treatment)+
        annotate("text", x = 45, y = 0.25, label = "Nutrients", size=3)+
        annotate("text", x = 99, y = 0.25, label = "Heat", size=3)
Figure3c

YII GLMs

All spp pooled (does not consider domminant symbiont)

# All spp together
      LME1<-lmer(YII ~ Treatment * DaysF * Spp + (1|Genotype),
                REML=TRUE, data=YII.data, na.action=na.omit)
       lmerTest::step (LME1)
## Backward reduced random-effect table:
## 
##                Eliminated npar logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                     154 8708.5 -17109                         
## (1 | Genotype)          0  153 8299.7 -16293 817.54  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                     Eliminated Sum Sq  Mean Sq NumDF  DenDF F value
## Treatment:DaysF:Spp          0 1.2893 0.020466    63 4851.1   14.11
##                        Pr(>F)    
## Treatment:DaysF:Spp < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment * DaysF * Spp + (1 | Genotype)
       drop1(LME1, test = "Chisq")
       #summary(LME1)
       anova(LME1)
       ranova(LME1)
      # EMMs
      All.YII.emm<-emmeans(LME1, ~Treatment * DaysF* Spp)
      emmip(LME1, ~DaysF|Treatment|Spp, CIs = TRUE) + theme_bw() + facet_grid(Spp~Treatment)

  # Spp responded differently, do separate analysis for each one

A.cer model

Nutrient treatment

Subset Acervicornis data

YII.Acer<-subset(YII.data, Spp=="Ac")
YII.Acer$Nutrients<-"Nutrients"
YII.Acer$Nutrients[YII.Acer$Treatment=="A"]<-"Ambient"

Find best model

# 1. Find the best model
YII.Acer$DaysF<-as.factor(YII.Acer$Days)

LME_Acer<-lmerTest::lmer(YII ~ Treatment * DaysF + 
                             (1|Genotype) + (1|Replicate) +  (1|Fragment), 
                              data=YII.Acer, na.action=na.omit)
      #summary(LME_Acer)
      
      Step.LME_Acer<-step (LME_Acer) # Replicate is not significant
      anova(LME_Acer)
      ranova(LME_Acer)# Replicate is not significant
      # Drop (1|Replicate)
      final_fm <- get_model(Step.LME_Acer)
      #summary(final_fm)
      
  LME_Acer1<-lmerTest::lmer(YII ~ Treatment * DaysF + 
                             (1|Genotype/Fragment), 
                              data=subset(YII.data, Spp=="Ac"), na.action=na.omit)
      
  ranova(LME_Acer1)
  LME_Acer2<-lmerTest::lmer(YII ~ Treatment * DaysF + 
                             (1|Genotype) +  (1|Fragment), 
                              data=subset(YII.data, Spp=="Ac"), na.action=na.omit)
  ranova(LME_Acer2)    
  anova(LME_Acer1, LME_Acer2) # LME_Acer1 and LME_Acer2 are the same 
#2. Extract EMMs
      Acer.YII.emm<-emmeans(LME_Acer1, ~Treatment | DaysF)
      contrast(Acer.YII.emm, "tukey")
## DaysF = 1:
##  contrast estimate      SE   df t.ratio p.value
##  A - N     0.01123 0.00566 1023  1.985  0.1164 
##  A - N+P   0.00684 0.00569 1022  1.202  0.4524 
##  N - N+P  -0.00439 0.00562 1023 -0.781  0.7147 
## 
## DaysF = 8:
##  contrast estimate      SE   df t.ratio p.value
##  A - N    -0.01626 0.00566 1023 -2.874  0.0115 
##  A - N+P  -0.01497 0.00569 1022 -2.630  0.0235 
##  N - N+P   0.00129 0.00562 1023  0.229  0.9715 
## 
## DaysF = 14:
##  contrast estimate      SE   df t.ratio p.value
##  A - N     0.00304 0.00566 1023  0.537  0.8533 
##  A - N+P  -0.00902 0.00569 1022 -1.584  0.2528 
##  N - N+P  -0.01206 0.00562 1023 -2.145  0.0815 
## 
## DaysF = 21:
##  contrast estimate      SE   df t.ratio p.value
##  A - N    -0.03309 0.00566 1023 -5.849  <.0001 
##  A - N+P  -0.02066 0.00569 1022 -3.629  0.0009 
##  N - N+P   0.01243 0.00562 1023  2.211  0.0698 
## 
## DaysF = 28:
##  contrast estimate      SE   df t.ratio p.value
##  A - N    -0.03049 0.00566 1023 -5.389  <.0001 
##  A - N+P  -0.03612 0.00569 1022 -6.344  <.0001 
##  N - N+P  -0.00562 0.00562 1023 -1.000  0.5768 
## 
## DaysF = 49:
##  contrast estimate      SE   df t.ratio p.value
##  A - N    -0.05631 0.00566 1023 -9.952  <.0001 
##  A - N+P  -0.05586 0.00569 1022 -9.811  <.0001 
##  N - N+P   0.00045 0.00562 1023  0.080  0.9965 
## 
## DaysF = 65:
##  contrast estimate      SE   df t.ratio p.value
##  A - N    -0.03088 0.00566 1023 -5.457  <.0001 
##  A - N+P  -0.04943 0.00573 1031 -8.633  <.0001 
##  N - N+P  -0.01855 0.00565 1032 -3.281  0.0031 
## 
## DaysF = 71:
##  contrast estimate      SE   df t.ratio p.value
##  A - N    -0.04150 0.00575 1050 -7.212  <.0001 
##  A - N+P  -0.02138 0.00573 1031 -3.734  0.0006 
##  N - N+P   0.02012 0.00575 1059  3.500  0.0014 
## 
## DaysF = 76:
##  contrast estimate      SE   df t.ratio p.value
##  A - N    -0.03999 0.00579 1061 -6.909  <.0001 
##  A - N+P  -0.03217 0.00573 1031 -5.619  <.0001 
##  N - N+P   0.00782 0.00578 1070  1.353  0.3664 
## 
## DaysF = 84:
##  contrast estimate      SE   df t.ratio p.value
##  A - N    -0.02435 0.00644 1226 -3.784  0.0005 
##  A - N+P  -0.03444 0.00634 1194 -5.436  <.0001 
##  N - N+P  -0.01009 0.00652 1252 -1.547  0.2692 
## 
## DaysF = 89:
##  contrast estimate      SE   df t.ratio p.value
##  A - N    -0.02714 0.00649 1238 -4.179  0.0001 
##  A - N+P  -0.02383 0.00634 1194 -3.760  0.0005 
##  N - N+P   0.00331 0.00658 1263  0.503  0.8698 
## 
## DaysF = 92:
##  contrast estimate      SE   df t.ratio p.value
##  A - N     0.01667 0.00676 1295  2.466  0.0367 
##  A - N+P  -0.01762 0.00649 1230 -2.713  0.0186 
##  N - N+P  -0.03429 0.00699 1340 -4.906  <.0001 
## 
## DaysF = 96:
##  contrast estimate      SE   df t.ratio p.value
##  A - N     0.05921 0.00712 1360  8.311  <.0001 
##  A - N+P   0.02681 0.00712 1352  3.764  0.0005 
##  N - N+P  -0.03240 0.00790 1452 -4.103  0.0001 
## 
## DaysF = 99:
##  contrast estimate      SE   df t.ratio p.value
##  A - N     0.12685 0.00736 1393 17.243  <.0001 
##  A - N+P   0.04875 0.00736 1388  6.627  <.0001 
##  N - N+P  -0.07810 0.00831 1485 -9.397  <.0001 
## 
## DaysF = 103:
##  contrast estimate      SE   df t.ratio p.value
##  A - N     0.12243 0.00781 1447 15.678  <.0001 
##  A - N+P   0.10504 0.00749 1406 14.022  <.0001 
##  N - N+P  -0.01739 0.00883 1516 -1.970  0.1200 
## 
## DaysF = 106:
##  contrast estimate      SE   df t.ratio p.value
##  A - N     0.18950 0.00877 1513 21.601  <.0001 
##  A - N+P   0.14946 0.00821 1478 18.205  <.0001 
##  N - N+P  -0.04004 0.01026 1555 -3.904  0.0003 
## 
## DaysF = 110:
##  contrast estimate      SE   df t.ratio p.value
##  A - N      nonEst      NA   NA     NA      NA 
##  A - N+P   0.13704 0.01142 1565 12.001  <.0001 
##  N - N+P    nonEst      NA   NA     NA      NA 
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
      Acer.YII.emm<-emmeans(LME_Acer1, ~Treatment * DaysF)
      
      # Effect plot options
      emmip(LME_Acer, ~DaysF|Treatment, CIs = TRUE) + theme_bw() # interaction plot of predictions

      Acer.YII_groups<-cld(Acer.YII.emm, by=NULL) # compact-letter display
      Acer.YII_groups<-Acer.YII_groups[order(Acer.YII_groups$Treatment, Acer.YII_groups$Day),]
      Acer.YII_groups
      #write.csv(Acer.YII_groups, "Outputs/Multicomp_AcerYII.csv", row.names = F)

Ofav and Ssid

Prepare data sets:

  • Time as discrete
#Time as a factor, not as int
      str(YII.data)
## 'data.frame':    5017 obs. of  17 variables:
##  $ Sample          : Factor w/ 6018 levels "Ac_101_T0","Ac_101_T1",..: 13 14 15 16 17 3 4 5 6 7 ...
##  $ Date            : Date, format: "2017-11-16" "2017-11-23" ...
##  $ Spp             : Factor w/ 3 levels "Ac","Of","Ss": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Fragment        : Factor w/ 354 levels "Ac_101","Ac_102",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Treatment       : Factor w/ 3 levels "A","N","N+P": 3 3 3 3 3 3 3 3 3 3 ...
##  $ Replicate       : Factor w/ 2 levels "R1","R2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ YII             : num  0.589 0.595 0.606 0.605 0.606 0.625 0.593 0.567 0.593 0.573 ...
##  $ Genotype        : Factor w/ 17 levels "G_48","G_62",..: 6 6 6 6 6 6 6 6 6 6 ...
##  $ Days            : num  1 8 14 21 28 49 65 71 76 84 ...
##  $ Time_Point      : Factor w/ 25 levels "T0","T1","T10",..: 21 22 23 24 25 3 4 5 6 7 ...
##  $ Phase           : Factor w/ 5 levels "Baseline","Heat",..: 1 3 3 3 3 3 3 3 3 4 ...
##  $ TotalSH         : num  0.1262 NA NA NA 0.0401 ...
##  $ logSH           : num  -0.899 NA NA NA -1.397 ...
##  $ D.Prp           : num  0 0 0 0 0 0 0 0 0 0 ...
##  $ Community       : Factor w/ 5 levels "A","B","C3","C1",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ InitialCommunity: Factor w/ 5 levels "A","B","C1","C3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ DaysF           : Factor w/ 26 levels "-112","-77","-36",..: 6 7 8 9 10 11 12 13 14 15 ...
      YII.data$DaysF<-as.factor(YII.data$Days)
  • Separate coral species

  • Separate phases:
    • BaseLine
    • C (nutrients only)
    • H (Heat challenge)

Ofav:

YII.Ofav<-subset(YII.data, Spp=="Of")

YII.Ofav.0<-subset(YII.Ofav, Days<2)
YII.Ofav.C<-subset(YII.Ofav, Days<77)
YII.Ofav.C.1<-subset(YII.Ofav.C, Days>2)
YII.Ofav.H<-subset(YII.Ofav, Days>75)

Ssid

YII.Ssid<-subset(YII.data, Spp=="Ss")

YII.Ssid.0<-subset(YII.Ssid, Days<2)
YII.Ssid.C<-subset(YII.Ssid, Days<77)
YII.Ssid.C.1<-subset(YII.Ssid.C, Days>2)
YII.Ssid.H<-subset(YII.Ssid, Days>75)

Baseline: effect of different Symbiodiniaceae taxa before treatments

Ofav (day 0)

LME_Ofav0<-lmer(YII ~ Treatment + InitialCommunity + (1|Replicate), 
                     data=YII.Ofav.0)
     #step (LME_Ofav0) 
     anova(LME_Ofav0) # Treatment not significant        
     ranova(LME_Ofav0) #  Replicate is not significant         
LME_Ofav0<-lm(YII ~ InitialCommunity, data=YII.Ofav.0)
     step (LME_Ofav0) 
## Start:  AIC=-508.95
## YII ~ InitialCommunity
## 
##                    Df Sum of Sq      RSS     AIC
## <none>                          0.057983 -508.95
## - InitialCommunity  1  0.014183 0.072166 -495.19
## 
## Call:
## lm(formula = YII ~ InitialCommunity, data = YII.Ofav.0)
## 
## Coefficients:
##       (Intercept)  InitialCommunityD  
##           0.53931           -0.03376
     anova(LME_Ofav0)         
# EMMs
    Ofav.YII.emm0<-emmeans(LME_Ofav0, ~ InitialCommunity)
      # contrast(Ssid.YII.emm, "tukey")
      
      # Tukey comparison (do not trust CI to compare EMMs)
      plot(emmeans(LME_Ofav0, ~InitialCommunity), comparisons = TRUE) +
        coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) +
        theme_bw()

      Ofav.YII_groups0<-cld(Ofav.YII.emm0, by=NULL) # compact-letter display
      Ofav.YII_groups0
      # write.csv(Ofav.YII_groups0, "Outputs/Multicomp_OfavYII0.csv", row.names = F)

Ssid (day 0)

LME_Ssid0<-lmer(YII ~ Treatment + InitialCommunity + (1|Replicate), 
                     data=YII.Ssid.0)
     step (LME_Ssid0) #  Treatemnt and replicate are significant :/
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                        7 369.23 -724.46                         
## (1 | Replicate)          0    6 338.52 -665.03 61.431  1  4.586e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                  Eliminated   Sum Sq  Mean Sq NumDF DenDF F value
## Treatment                 0 0.019079 0.009540     2   156  20.809
## InitialCommunity          0 0.087511 0.043756     2   156  95.445
##                     Pr(>F)    
## Treatment        9.754e-09 ***
## InitialCommunity < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + InitialCommunity + (1 | Replicate)
     anova(LME_Ssid0)         
     ranova(LME_Ssid0)         
# EMMs
      Ssid.YII.emm0<-emmeans(LME_Ssid0, ~ InitialCommunity)
      # contrast(Ssid.YII.emm, "tukey")
      
      # Tukey comparison (do not trust CI to compare EMMs)
      plot(emmeans(LME_Ssid0, ~InitialCommunity), comparisons = TRUE) +
        coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) +
        theme_bw()

      Ssid.YII_groups0<-cld(Ssid.YII.emm0, by=NULL) # compact-letter display
      Ssid.YII_groups0
      # write.csv(Ssid.YII_groups0, "Outputs/Multicomp_SsidYII0.csv", row.names = F)

C: Effect of nutrient treatments at control temperature

Ofav

  • C Days pooled (days 8-76)
LME_Ofav1.1<-lmer(YII ~ Treatment * InitialCommunity + (1|Fragment), 
                     data=YII.Ofav.C.1)
     step (LME_Ofav1.1) #  Replicate is significant
## Backward reduced random-effect table:
## 
##                Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                       8 1046.0 -2075.9                         
## (1 | Fragment)          0    7 1036.4 -2058.8 19.116  1   1.23e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                            Eliminated   Sum Sq  Mean Sq NumDF DenDF
## Treatment:InitialCommunity          1 0.002297 0.001148     2    66
## Treatment                           0 0.013198 0.006599     2    68
## InitialCommunity                    0 0.173819 0.173819     1    68
##                             F value    Pr(>F)    
## Treatment:InitialCommunity   0.8780  0.420421    
## Treatment                    5.0455  0.009055 ** 
## InitialCommunity           132.8935 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + InitialCommunity + (1 | Fragment)
     anova(LME_Ofav1.1)         
     ranova(LME_Ofav1.1)         
# EMMs
    Ofav.YII.emm1.1<-emmeans(LME_Ofav1.1, ~ InitialCommunity|Treatment)
    # contrast(Ssid.YII.emm, "tukey")
      
    # Tukey comparison (do not trust CI to compare EMMs)
      plot(emmeans(LME_Ofav1.1, ~Treatment|InitialCommunity), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_grid(~InitialCommunity) +
        theme_bw()

      Ofav.YII.emm1.1<-cld(Ofav.YII.emm1.1, by=NULL) # compact-letter display
      Ofav.YII.emm1.1
      #write.csv(Ofav.YII.emm1.1, "Outputs/Multicomp_OfavYII1.csv", row.names = F)
  • Model for “Days” as continous variable (days 1-76)
# 1. Ofav- Days 
      LME_Ofav<-lmer(YII ~ Treatment * Days * InitialCommunity +
                       (1|Genotype) + (1|Fragment) +(1|Replicate), 
                        data=YII.Ofav.C, na.action=na.omit)
      lmerTest::step (LME_Ofav) # Replicate is not significant
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                       16 1243.0 -2454.1                         
## (1 | Replicate)          1   15 1243.0 -2456.1  0.000  1          1    
## (1 | Genotype)           0   14 1233.7 -2439.4 18.642  1  1.577e-05 ***
## (1 | Fragment)           0   14 1233.5 -2438.9 19.149  1  1.209e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                 Eliminated    Sum Sq   Mean Sq NumDF DenDF
## Treatment:Days:InitialCommunity          0 0.0095086 0.0047543     2   570
##                                 F value   Pr(>F)   
## Treatment:Days:InitialCommunity   5.128 0.006205 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + Days + InitialCommunity + (1 | Genotype) + 
##     (1 | Fragment) + Treatment:Days + Treatment:InitialCommunity + 
##     Days:InitialCommunity + Treatment:Days:InitialCommunity
      ranova(LME_Ofav)
      LME_Ofav1<-lmer(YII ~ Treatment * Days * InitialCommunity + 
                        (1 | Genotype/Fragment),
                      data=YII.Ofav.C, na.action=na.omit)
      lmerTest::step (LME_Ofav1)
## Backward reduced random-effect table:
## 
##                         Eliminated npar logLik     AIC    LRT Df
## <none>                               15 1243.0 -2456.1          
## (1 | Fragment:Genotype)          0   14 1233.5 -2438.9 19.149  1
## (1 | Genotype)                   0   14 1233.7 -2439.4 18.642  1
##                         Pr(>Chisq)    
## <none>                                
## (1 | Fragment:Genotype)  1.209e-05 ***
## (1 | Genotype)           1.577e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                 Eliminated    Sum Sq   Mean Sq NumDF DenDF
## Treatment:Days:InitialCommunity          0 0.0095086 0.0047543     2   570
##                                 F value   Pr(>F)   
## Treatment:Days:InitialCommunity   5.128 0.006205 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment * Days * InitialCommunity + (1 | Genotype/Fragment)
# 2. Predict values:
    pred_Ofav1 <- predict(LME_Ofav1,re.form = NA)

  #3. Bootstrap CI:
    OF.boot1 <- bootMer(LME_Ofav1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(OF.boot1$t, 2, sd)
    CI.lo_1 <- pred_Ofav1 - std.err*1.96
    CI.hi_1 <- pred_Ofav1 + std.err*1.96

  #Plot
  Model_of_1b_plot<- ggplot(
    YII.Ofav.C, aes(x = Days, y = YII, colour = Treatment)) +
    geom_line(aes(y = pred_Ofav1),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    scale_y_continuous(name=expression(~italic("Fv / Fm")),
                      limits = c(0.35, 0.61), 
                      breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
    scale_x_continuous("Days in the experiment", limits = c(0, 78),
                     breaks = seq(0, 76, by=7), expand = c(0,0))+
    
     stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                linetype=1, alpha=0.5) + 

    # stat_summary(fun.y=mean, geom="point", size =1,
    #             position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_of_1b_plot + facet_grid(~InitialCommunity)

  • Model for “Days” as factor
# 1. Ofav- Days 
      LME_Ofav<-lmer(YII ~ Treatment * DaysF * InitialCommunity +
                       (1|Genotype) + (1|Fragment) +(1|Replicate), 
                        data=YII.Ofav.C, na.action=na.omit)
      lmerTest::step (LME_Ofav) # Replicate is not significant
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                       58 1311.5 -2507.1                         
## (1 | Replicate)          1   57 1311.5 -2509.1  0.000  1          1    
## (1 | Genotype)           0   56 1302.2 -2492.5 18.642  1  1.577e-05 ***
## (1 | Fragment)           0   56 1273.0 -2434.0 77.139  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                  Eliminated   Sum Sq   Mean Sq NumDF DenDF
## Treatment:DaysF:InitialCommunity          0 0.021644 0.0013528    16   528
##                                  F value    Pr(>F)    
## Treatment:DaysF:InitialCommunity  2.7855 0.0002475 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + DaysF + InitialCommunity + (1 | Genotype) + 
##     (1 | Fragment) + Treatment:DaysF + Treatment:InitialCommunity + 
##     DaysF:InitialCommunity + Treatment:DaysF:InitialCommunity
      ranova(LME_Ofav)
      LME_Ofav1.2<-lmer(YII ~ Treatment * DaysF * InitialCommunity + 
                        (1|Fragment),
                      data=YII.Ofav.C, na.action=na.omit)
      lmerTest::step (LME_Ofav1.2)
## Backward reduced random-effect table:
## 
##                Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                      56 1302.2 -2492.5                         
## (1 | Fragment)          0   55 1230.0 -2349.9 144.55  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                  Eliminated   Sum Sq   Mean Sq NumDF DenDF
## Treatment:DaysF:InitialCommunity          0 0.021644 0.0013528    16   528
##                                  F value    Pr(>F)    
## Treatment:DaysF:InitialCommunity  2.7856 0.0002475 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment * DaysF * InitialCommunity + (1 | Fragment)
  # 2. Predict values:
    pred_Ofav1.2 <- predict(LME_Ofav1.2,re.form = NA)

  #3. Bootstrap CI:
    OF.boot1.2 <- bootMer(LME_Ofav1.2, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(OF.boot1.2$t, 2, sd)
    CI.lo_1 <- pred_Ofav1.2 - std.err*1.96
    CI.hi_1 <- pred_Ofav1.2 + std.err*1.96

  #Plot
  Model_of_1.2_plot<- ggplot(
    YII.Ofav.C, aes(x = Days, y = YII, colour = Treatment)) +
    geom_line(aes(y = pred_Ofav1.2),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    scale_y_continuous(name=expression(~italic("Fv / Fm")),
                       limits = c(0.35, 0.61), 
                       breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
    scale_x_continuous("Days in the experiment", limits = c(0, 78),
                     breaks = seq(0, 76, by=7), expand = c(0,0))+
    
     stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
      stat_summary(fun.y=mean, geom="point", size =1,
                   position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_of_1.2_plot + facet_grid(~InitialCommunity)

  # EMMs
    Ofav.YII.emmC<-emmeans(LME_Ofav1.2, ~Treatment * DaysF * InitialCommunity)
    #contrast(Ofav.YII.emm, "tukey")
      
      
      # Ofav.YII_groupsC<-cld(Ofav.YII.emmC, by=NULL) # compact-letter display
      # Ofav.YII_groupsC<-Ofav.YII_groupsC[order(
      #   Ofav.YII_groupsC$Days,
      #   Ofav.YII_groupsC$Treatment,
      #   Ofav.YII_groupsC$InitialCommunity),]
      # Ofav.YII_groupsC
      # write.csv(Ofav.YII_groupsC, "Outputs/Multicomp_OfavYIIC.csv", row.names = F)

Ssid

  • Days pooled (8-76 Days)
LME_Ssid1.1<-lmer(YII ~ Treatment * InitialCommunity + (1|Fragment), 
                     data=YII.Ssid.C.1)
     step (LME_Ssid1.1) #  Replicate is significant
## Backward reduced random-effect table:
## 
##                Eliminated npar logLik     AIC   LRT Df Pr(>Chisq)    
## <none>                      11 2224.0 -4425.9                        
## (1 | Fragment)          0   10 2211.5 -4403.1 24.84  1  6.229e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                            Eliminated   Sum Sq   Mean Sq NumDF  DenDF
## Treatment:InitialCommunity          0 0.039046 0.0097615     4 153.16
##                            F value   Pr(>F)    
## Treatment:InitialCommunity  5.9109 0.000189 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment * InitialCommunity + (1 | Fragment)
     anova(LME_Ssid1.1)         
     ranova(LME_Ssid1.1)         
# EMMs
    Ssid.YII.emm1.1<-emmeans(LME_Ssid1.1, ~ InitialCommunity|Treatment)
    # contrast(Ssid.YII.emm, "tukey")
      
    # Tukey comparison (do not trust CI to compare EMMs)
      plot(emmeans(LME_Ssid1.1, ~Treatment|InitialCommunity), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_grid(~InitialCommunity)+
        theme_bw()

      Ssid.YII.Groups1.1<-cld(Ssid.YII.emm1.1, by=NULL) # compact-letter display
      Ssid.YII.Groups1.1<- Ssid.YII.Groups1.1[order(
        Ssid.YII.Groups1.1$InitialCommunity, 
        Ssid.YII.Groups1.1$Treatment),]
      Ssid.YII.Groups1.1
     #write.csv( Ssid.YII.Groups1.1, "Outputs/Multicomp_SsidYII1.1.csv", row.names = F)
  • Model for “Days” as continous variable (days 1-76)
LME_Ssid<-lmer(YII ~ Treatment * Days * InitialCommunity +
                       (1|Fragment) +(1|Replicate), 
                     data=subset(YII.Ssid.C))
     step (LME_Ssid) #  Replicate is not significant
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                       21 2486.7 -4931.3                         
## (1 | Replicate)          1   20 2486.7 -4933.3  0.000  1          1    
## (1 | Fragment)           0   19 2466.3 -4894.6 40.775  1  1.708e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                 Eliminated   Sum Sq   Mean Sq NumDF
## Treatment:Days:InitialCommunity          1 0.004115 0.0010287     4
## Treatment:Days                           0 0.028806 0.0144029     2
## Treatment:InitialCommunity               0 0.033666 0.0084166     4
## Days:InitialCommunity                    0 0.015350 0.0076750     2
##                                   DenDF F value    Pr(>F)    
## Treatment:Days:InitialCommunity 1283.90  0.6706 0.6124001    
## Treatment:Days                  1288.52  9.3992 8.861e-05 ***
## Treatment:InitialCommunity       153.14  5.4926 0.0003694 ***
## Days:InitialCommunity           1287.90  5.0087 0.0068106 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + Days + InitialCommunity + (1 | Fragment) + 
##     Treatment:Days + Treatment:InitialCommunity + Days:InitialCommunity
     anova(LME_Ssid)         
     ranova(LME_Ssid)         
LME_Ssid1<-lmer(YII ~ Treatment * Days * InitialCommunity +
                       (1|Fragment),
                     data=YII.Ssid.C, na.action=na.omit)


# EMMs
      Ssid.YII.emm1<-emmeans(LME_Ssid1, ~Treatment * Days * InitialCommunity)
      # contrast(Ssid.YII.emm, "tukey")
      
      # Effect plot options
      emmip(LME_Ssid1, ~InitialCommunity|Treatment, CIs = TRUE) + theme_bw() # interaction plot of predictions

# 2. Predict values:
    predSs1 <- predict(LME_Ssid1,re.form = NA)

# 3. Bootstrap CI:
    boot1Ss <- bootMer(LME_Ssid1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(boot1Ss$t, 2, sd)
    CI.lo_1 <- predSs1 - std.err*1.96
    CI.hi_1 <- predSs1 + std.err*1.96
# 4. Plot
  Model1Ss_plot<- ggplot(
    YII.Ssid.C, aes(x = Days, y = YII, colour = Treatment)) +
    geom_line(aes(y = predSs1),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.5) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    scale_y_continuous(name=expression(~italic("Fv / Fm")),
                       limits = c(0.3, 0.6), 
                       breaks = seq(0.3, 0.6, by=0.1), expand = c(0,0))+
    scale_x_continuous("Days in the experiment (C phase)", limits = c(0, 78),
                     breaks = seq(0, 76, by=7), expand = c(0,0))+
    
    stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                   linetype=1, alpha=0.5) + 
    #stat_summary(fun.y=mean, geom="point", size =1,
    #               position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw # + Fill.colour
  
  Model1Ss_plot + facet_grid(~InitialCommunity)

  • Model for “Days” as factor
# 1. Ssid- Days 
      LME_Ssid<-lmer(YII ~ Treatment * DaysF * InitialCommunity +
                       (1|Genotype) + (1|Fragment) +(1|Replicate), 
                        data=YII.Ssid.C, na.action=na.omit)
      lmerTest::step (LME_Ssid) # Replicate is not significant
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC     LRT Df Pr(>Chisq)    
## <none>                       85 3033.7 -5897.4                          
## (1 | Replicate)          1   84 3033.3 -5898.6   0.841  1      0.359    
## (1 | Genotype)           0   83 2995.3 -5824.7  75.896  1     <2e-16 ***
## (1 | Fragment)           0   83 2972.7 -5779.4 121.188  1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                  Eliminated   Sum Sq   Mean Sq NumDF
## Treatment:DaysF:InitialCommunity          1 0.019444 0.0006076    32
## Treatment:DaysF                           0 0.134201 0.0083875    16
## Treatment:InitialCommunity                0 0.020396 0.0050991     4
## DaysF:InitialCommunity                    0 0.044362 0.0027727    16
##                                    DenDF F value    Pr(>F)    
## Treatment:DaysF:InitialCommunity 1220.34  1.1718    0.2354    
## Treatment:DaysF                  1252.49 16.1060 < 2.2e-16 ***
## Treatment:InitialCommunity        147.69  9.7915 4.712e-07 ***
## DaysF:InitialCommunity           1252.33  5.3241 5.033e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + DaysF + InitialCommunity + (1 | Genotype) + 
##     (1 | Fragment) + Treatment:DaysF + Treatment:InitialCommunity + 
##     DaysF:InitialCommunity
      ranova(LME_Ssid)
      LME_Ssid1.2<-lmer(YII ~ Treatment * DaysF * InitialCommunity + 
                        (1|Fragment),
                      data=YII.Ssid.C, na.action=na.omit)
      lmerTest::step (LME_Ssid1.2)
## Backward reduced random-effect table:
## 
##                Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                      83 2995.3 -5824.7                         
## (1 | Fragment)          0   82 2829.0 -5494.0 332.65  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                  Eliminated   Sum Sq   Mean Sq NumDF
## Treatment:DaysF:InitialCommunity          1 0.019389 0.0006059    32
## Treatment:DaysF                           0 0.134162 0.0083851    16
## Treatment:InitialCommunity                0 0.011521 0.0028803     4
## DaysF:InitialCommunity                    0 0.044336 0.0027710    16
##                                    DenDF F value    Pr(>F)    
## Treatment:DaysF:InitialCommunity 1220.22  1.1685 0.2390861    
## Treatment:DaysF                  1252.32 16.1009 < 2.2e-16 ***
## Treatment:InitialCommunity        153.05  5.5307 0.0003475 ***
## DaysF:InitialCommunity           1252.23  5.3208  5.14e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + DaysF + InitialCommunity + (1 | Fragment) + 
##     Treatment:DaysF + Treatment:InitialCommunity + DaysF:InitialCommunity
# 2. Predict values:
    pred_Ssid1.2 <- predict(LME_Ssid1.2,re.form = NA)

  #3. Bootstrap CI:
    Ss.boot1.2 <- bootMer(LME_Ssid1.2, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(Ss.boot1.2$t, 2, sd)
    CI.lo_1 <- pred_Ssid1.2 - std.err*1.96
    CI.hi_1 <- pred_Ssid1.2 + std.err*1.96
# 4. Plot
  Model_Ss_1.2_plot<- ggplot(
    YII.Ssid.C, aes(x = Days, y = YII, colour = Treatment)) +
    geom_line(aes(y = pred_Ssid1.2),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    scale_y_continuous(name=expression(~italic("Fv / Fm")),
                       limits = c(0.35, 0.6), 
                       breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
    scale_x_continuous("Days in the experiment", limits = c(0, 78),
                     breaks = seq(0, 76, by=7), expand = c(0,0))+
    
    stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="point", size =1,
                   position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_Ss_1.2_plot + facet_grid(~InitialCommunity)

# 5. EMMs
    Ssid.YII.emmC<-emmeans(LME_Ssid1.2, ~Treatment * DaysF * InitialCommunity)
      #contrast(Ssid.YII.emm, "tukey")
      
      
       Ssid.YII_groupsC<-cld(Ssid.YII.emmC, by=NULL) # compact-letter display
       Ssid.YII_groupsC<-Ssid.YII_groupsC[order(
         Ssid.YII_groupsC$DaysF,
         Ssid.YII_groupsC$InitialCommunity,
         Ssid.YII_groupsC$Treatment),]
       Ssid.YII_groupsC
       #write.csv(Ssid.YII_groupsC, "Outputs/Multicomp_SsidYIIC.csv", row.names = F)

H: Effect of pre-exposure to nutrient treatments during heat challenge

Ofav

  • Model for “Days” as continous variable (days 84-110)
# 1. Ofav- Days 
      LME_Ofav.H<-lmer(YII ~ Treatment * Days * InitialCommunity +
                       (1|Genotype) + (1|Fragment) +(1|Replicate), 
                        data=YII.Ofav.H, na.action=na.omit)
      lmerTest::step (LME_Ofav.H) # Replicate is not significant
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)   
## <none>                       16 666.73 -1301.5                        
## (1 | Replicate)          1   15 666.73 -1303.5 0.0000  1   1.000000   
## (1 | Genotype)           0   14 662.44 -1296.9 8.5781  1   0.003402 **
## (1 | Fragment)           0   14 663.16 -1298.3 7.1280  1   0.007589 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                 Eliminated    Sum Sq   Mean Sq NumDF
## Treatment:Days:InitialCommunity          1 0.0022744 0.0011372     2
## Treatment:InitialCommunity               2 0.0016208 0.0008104     2
## Treatment:Days                           0 0.0189804 0.0094902     2
## Days:InitialCommunity                    0 0.0088774 0.0088774     1
##                                  DenDF F value   Pr(>F)   
## Treatment:Days:InitialCommunity 390.99  0.7378 0.478849   
## Treatment:InitialCommunity       89.11  0.5261 0.592744   
## Treatment:Days                  390.47  6.1692 0.002302 **
## Days:InitialCommunity           393.25  5.7709 0.016757 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + Days + InitialCommunity + (1 | Genotype) + 
##     (1 | Fragment) + Treatment:Days + Days:InitialCommunity
      ranova(LME_Ofav.H)
      LME_Ofav.H1<-lmer(YII ~ Treatment * Days + 
                        (1 | Genotype/Fragment),
                      data=YII.Ofav.H, na.action=na.omit)
      lmerTest::step (LME_Ofav.H1)
## Backward reduced random-effect table:
## 
##                         Eliminated npar logLik     AIC     LRT Df
## <none>                                9 687.23 -1356.5           
## (1 | Fragment:Genotype)          0    8 683.59 -1351.2  7.2867  1
## (1 | Genotype)                   0    8 674.23 -1332.5 26.0040  1
##                         Pr(>Chisq)    
## <none>                                
## (1 | Fragment:Genotype)   0.006947 ** 
## (1 | Genotype)           3.407e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                Eliminated   Sum Sq   Mean Sq NumDF DenDF F value  Pr(>F)
## Treatment:Days          0 0.015625 0.0078127     2   390  5.0081 0.00712
##                  
## Treatment:Days **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment * Days + (1 | Genotype/Fragment)
# 2. Predict values:
    pred_Ofav.H <- predict(LME_Ofav.H1,re.form = NA)

#3. Bootstrap CI:
    OF.bootH <- bootMer(LME_Ofav.H1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(OF.bootH$t, 2, sd)
    CI.lo_H <- pred_Ofav.H - std.err*1.96
    CI.hi_H <- pred_Ofav.H + std.err*1.96
 # 4 .Plot
  Model_of_H_plot<- ggplot(
    YII.Ofav.H, aes(x = Days, y = YII, colour = Treatment)) +
    geom_line(aes(y = pred_Ofav.H),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_H, ymax = CI.hi_H),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    scale_y_continuous(name=expression(~italic("Fv / Fm")),
                       limits = c(0.25, 0.58), 
                       breaks = seq(0.3, 0.5, by=0.1), expand = c(0,0))+
    scale_x_continuous("Days in the experiment (H phase)", limits = c(80, 112),
                     breaks = seq(80, 110, by=7), expand = c(0,0))+
    
    stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                   linetype=1, alpha=1) + 
    stat_summary(fun.y=mean, geom="point", size =1,
                   position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_of_H_plot

  Model_of_H_plot+ facet_grid(InitialCommunity~.)

  • Model for “Days” as factor (daysF 84-110)
# 1. Ofav- Days 
    LME_Ofav.HF<-lmer(YII ~ Treatment * DaysF * Community +
                       (1|Genotype) + (1|Fragment) +(1|Replicate), 
                        data=YII.Ofav.H, na.action=na.omit)
      lmerTest::step (LME_Ofav.HF) # Replicate is not significant
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                       49 740.58 -1383.2                         
## (1 | Replicate)          1   48 740.58 -1385.2  0.000  1   1.000000    
## (1 | Genotype)           0   47 735.75 -1377.5  9.650  1   0.001894 ** 
## (1 | Fragment)           0   47 715.64 -1337.3 49.865  1  1.647e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated   Sum Sq   Mean Sq NumDF  DenDF
## Treatment:DaysF:Community          1 0.001743 0.0002179     8 330.04
## Treatment:Community                2 0.000134 0.0000672     2 110.06
## Treatment:DaysF                    0 0.033700 0.0021062    16 330.20
## DaysF:Community                    0 0.047109 0.0067299     7 336.98
##                           F value    Pr(>F)    
## Treatment:DaysF:Community  0.3417    0.9492    
## Treatment:Community        0.1071    0.8986    
## Treatment:DaysF            3.3619 1.663e-05 ***
## DaysF:Community           10.7421 3.153e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + DaysF + Community + (1 | Genotype) + (1 | Fragment) + 
##     Treatment:DaysF + DaysF:Community
      ranova(LME_Ofav.HF)
      LME_Ofav.HF.1<-lmer(YII ~ Treatment + DaysF + Community + 
                          (1 | Genotype) + (1 | Fragment) +
                          Treatment:DaysF + DaysF:Community,
                      data=YII.Ofav.H, na.action=na.omit)
      lmerTest::step (LME_Ofav.HF.1)
## Backward reduced random-effect table:
## 
##                Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                      38 767.61 -1459.2                         
## (1 | Genotype)          0   37 762.53 -1451.1 10.158  1   0.001437 ** 
## (1 | Fragment)          0   37 742.39 -1410.8 50.450  1  1.222e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                 Eliminated   Sum Sq   Mean Sq NumDF  DenDF F value
## Treatment:DaysF          0 0.033700 0.0021062    16 330.20  3.3619
## DaysF:Community          0 0.047109 0.0067299     7 336.98 10.7421
##                    Pr(>F)    
## Treatment:DaysF 1.663e-05 ***
## DaysF:Community 3.153e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + DaysF + Community + (1 | Genotype) + (1 | Fragment) + 
##     Treatment:DaysF + DaysF:Community
# 2. Predict values:
    pred_OfavHF <- predict(LME_Ofav.HF.1,re.form = NA)

# 3. Bootstrap CI:
    OF.bootHF <- bootMer(LME_Ofav.HF.1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(OF.bootHF$t, 2, sd)
    CI.lo_HF <- pred_OfavHF - std.err*1.96
    CI.hi_HF <- pred_OfavHF + std.err*1.96
# 4. Plot
  Model_of_HF_plot<- ggplot(
    YII.Ofav.H, aes(x = Days, y = YII, colour = Treatment)) +
    geom_line(aes(y = pred_OfavHF),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_HF, ymax = CI.hi_HF),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    scale_y_continuous(name=expression(~italic("Fv / Fm")), 
                       limits = c(0.25, 0.55), 
                       breaks = seq(0.3, 0.5, by=0.1), expand = c(0,0))+
    scale_x_continuous("Days in the experiment (H phase)", limits = c(80, 112),
                     breaks = seq(80, 110, by=7), expand = c(0,0))+
    
    stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="point", size =1,
                   position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_of_HF_plot

  Model_of_HF_plot + facet_grid(Community~.)

# 5. EMMs
      Ofav.YII.emmH<-emmeans(LME_Ofav.HF.1, ~ Treatment:DaysF + DaysF:Community)
      #contrast(Ofav.YII.emmH, "tukey")
      
      Ofav.YII_groupsH<-cld(Ofav.YII.emmH, by=NULL) # compact-letter display
      Ofav.YII_groupsH<-Ofav.YII_groupsH[order(
        Ofav.YII_groupsH$Days,
        Ofav.YII_groupsH$Treatment,
        Ofav.YII_groupsH$Community),]
      Ofav.YII_groupsH
      #write.csv(Ofav.YII_groupsH, "Outputs/Multicomp_OfavYII_H.csv", row.names = F)

Ssid

  • Model for “Days” as continous variable (days 84-110)
# 1. Ssid- Days 
      LME_Ssid.H<-lmer(YII ~ Treatment * Days * InitialCommunity +
                       (1|Genotype) + (1|Fragment) +(1|Replicate), 
                        data=YII.Ssid.H, na.action=na.omit)
      lmerTest::step (LME_Ssid.H) # Replicate is not significant
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)   
## <none>                       22 1692.9 -3341.9                        
## (1 | Replicate)          1   21 1692.9 -3343.9 0.0000  1   1.000000   
## (1 | Genotype)           0   20 1690.8 -3341.5 4.3147  1   0.037785 * 
## (1 | Fragment)           0   20 1688.9 -3337.9 7.9765  1   0.004739 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                 Eliminated   Sum Sq  Mean Sq NumDF DenDF
## Treatment:Days:InitialCommunity          0 0.043395 0.010849     4  1029
##                                 F value    Pr(>F)    
## Treatment:Days:InitialCommunity  4.6856 0.0009452 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + Days + InitialCommunity + (1 | Genotype) + 
##     (1 | Fragment) + Treatment:Days + Treatment:InitialCommunity + 
##     Days:InitialCommunity + Treatment:Days:InitialCommunity
      ranova(LME_Ssid.H)
      LME_Ssid.H1<-lmer(YII ~ Treatment * Days * InitialCommunity + 
                        (1 | Genotype/Fragment),
                      data=YII.Ssid.H, na.action=na.omit)
      lmerTest::step (LME_Ssid.H1)
## Backward reduced random-effect table:
## 
##                         Eliminated npar logLik     AIC    LRT Df
## <none>                               21 1692.9 -3343.9          
## (1 | Fragment:Genotype)          0   20 1688.9 -3337.9 7.9765  1
## (1 | Genotype)                   0   20 1690.8 -3341.5 4.3147  1
##                         Pr(>Chisq)   
## <none>                               
## (1 | Fragment:Genotype)   0.004739 **
## (1 | Genotype)            0.037785 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                 Eliminated   Sum Sq  Mean Sq NumDF DenDF
## Treatment:Days:InitialCommunity          0 0.043395 0.010849     4  1029
##                                 F value    Pr(>F)    
## Treatment:Days:InitialCommunity  4.6856 0.0009452 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment * Days * InitialCommunity + (1 | Genotype/Fragment)
      anova(LME_Ssid.H1)
# 2. Predict values:
    pred_Ssid.H <- predict(LME_Ssid.H1,re.form = NA)

#3. Bootstrap CI:
    Ss.bootH <- bootMer(LME_Ssid.H1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.errSH <- apply(Ss.bootH$t, 2, sd)
    CI.lo_SH <- pred_Ssid.H - std.errSH*1.96
    CI.hi_SH <- pred_Ssid.H + std.errSH*1.96
 # 4 .Plot
  Model_Ss_H_plot<- ggplot(
    YII.Ssid.H, aes(x = Days, y = YII, colour = Treatment)) +
    geom_line(aes(y = pred_Ssid.H),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_SH, ymax = CI.hi_SH),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    scale_y_continuous(name=expression(~italic("Fv / Fm")),
                       limits = c(0.1, 0.58), 
                       breaks = seq(0.1, 0.5, by=0.1), expand = c(0,0))+
    scale_x_continuous("Days in the experiment (H phase)", limits = c(80, 112),
                     breaks = seq(80, 110, by=7), expand = c(0,0))+
    
    stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                   linetype=1, alpha=0.5) + 
    #stat_summary(fun.y=mean, geom="point", size =1,
    #               position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_Ss_H_plot+ facet_grid(~InitialCommunity)

  • Model for “Days” as factor (daysF 84-110)
# 1. Ssid- Days 
      LME_Ssid.HF<-lmer(YII ~ Treatment * DaysF * InitialCommunity +
                       (1|Genotype) + (1|Fragment) +(1|Replicate), 
                        data=YII.Ssid.H, na.action=na.omit)
      lmerTest::step (LME_Ssid.HF) # Replicate is not significant
## Backward reduced random-effect table:
## 
##                 Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                       85 1885.7 -3601.5                         
## (1 | Replicate)          1   84 1885.7 -3603.5 -0.010  1     1.0000    
## (1 | Genotype)           0   83 1883.2 -3600.3  5.153  1     0.0232 *  
## (1 | Fragment)           0   83 1846.5 -3526.9 78.534  1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                  Eliminated   Sum Sq   Mean Sq NumDF
## Treatment:DaysF:InitialCommunity          0 0.081998 0.0025624    32
##                                   DenDF F value   Pr(>F)    
## Treatment:DaysF:InitialCommunity 929.23   2.275 8.27e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment + DaysF + InitialCommunity + (1 | Genotype) + 
##     (1 | Fragment) + Treatment:DaysF + Treatment:InitialCommunity + 
##     DaysF:InitialCommunity + Treatment:DaysF:InitialCommunity
      ranova(LME_Ssid.HF)
      LME_Ssid.HF.1<-lmer(YII ~ Treatment * DaysF * InitialCommunity + (1 | Fragment),
                      data=YII.Ssid.H, na.action=na.omit)
      lmerTest::step (LME_Ssid.HF.1)
## Backward reduced random-effect table:
## 
##                Eliminated npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>                      83 1883.2 -3600.3                         
## (1 | Fragment)          0   82 1835.2 -3506.5 95.807  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                                  Eliminated   Sum Sq   Mean Sq NumDF
## Treatment:DaysF:InitialCommunity          0 0.081916 0.0025599    32
##                                   DenDF F value    Pr(>F)    
## Treatment:DaysF:InitialCommunity 928.77  2.2693 8.706e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## YII ~ Treatment * DaysF * InitialCommunity + (1 | Fragment)
      anova(LME_Ssid.HF.1)
# 2. Predict values:
    pred_SsidHF <- predict(LME_Ssid.HF.1,re.form = NA)

# 3. Bootstrap CI:
    Ss.bootHF <- bootMer(LME_Ssid.HF.1, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(Ss.bootHF$t, 2, sd)
    CI.lo_HF.S <- pred_SsidHF - std.err*1.96
    CI.hi_HF.S <- pred_SsidHF + std.err*1.96
# 4. Plot
  Model_Ss_HF_plot<- ggplot(
    YII.Ssid.H, aes(x = Days, y = YII, colour = Treatment)) +
    geom_line(aes(y = pred_SsidHF),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_HF.S, ymax = CI.hi_HF.S),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    scale_y_continuous(name=expression(~italic("Fv / Fm")), 
                       limits = c(0.1, 0.55), 
                       breaks = seq(0.1, 0.5, by=0.1), expand = c(0,0))+
    scale_x_continuous("Days in the experiment (H phase)", limits = c(80, 112),
                     breaks = seq(80, 110, by=7), expand = c(0,0))+
    
    stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    #stat_summary(fun.y=mean, geom="point", size =1,
    #               position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_Ss_HF_plot + facet_grid(~InitialCommunity)

# 5. EMMs
    Ssid.YII.emmH<-emmeans(LME_Ssid.HF.1, ~ Treatment*DaysF*InitialCommunity)
      #contrast(Ssid.YII.emmH, "tukey")
      
      Ssid.YII_groupsH<-cld(Ssid.YII.emmH, by=NULL) # compact-letter display
      Ssid.YII_groupsH<-Ssid.YII_groupsH[order(
        Ssid.YII_groupsH$Days,
        Ssid.YII_groupsH$Treatment,
        Ssid.YII_groupsH$InitialCommunity),]
      Ssid.YII_groupsH
      #write.csv(Ssid.YII_groupsH, "Outputs/Multicomp_SsidYIIH.csv", row.names = F)

Packages used

# Creates bibliography 
#knitr::write_bib(c(.packages()), "packages.bib")

Arnold, Jeffrey B. 2019. Ggthemes: Extra Themes, Scales and Geoms for ’Ggplot2’. https://CRAN.R-project.org/package=ggthemes.

Bates, Douglas, and Martin Maechler. 2019. Matrix: Sparse and Dense Matrix Classes and Methods. https://CRAN.R-project.org/package=Matrix.

Bates, Douglas, Martin Maechler, Ben Bolker, and Steven Walker. 2019. Lme4: Linear Mixed-Effects Models Using ’Eigen’ and S4. https://CRAN.R-project.org/package=lme4.

Fox, John, Sanford Weisberg, and Brad Price. 2018. CarData: Companion to Applied Regression Data Sets. https://CRAN.R-project.org/package=carData.

Fox, John, Sanford Weisberg, Brad Price, Michael Friendly, and Jangman Hong. 2019. Effects: Effect Displays for Linear, Generalized Linear, and Other Models. https://CRAN.R-project.org/package=effects.

Genz, Alan, Frank Bretz, Tetsuhisa Miwa, Xuefei Mi, and Torsten Hothorn. 2019. Mvtnorm: Multivariate Normal and T Distributions. https://CRAN.R-project.org/package=mvtnorm.

Gohel, David, Hadley Wickham, Lionel Henry, and Jeroen Ooms. 2019. Gdtools: Utilities for Graphical Rendering. https://CRAN.R-project.org/package=gdtools.

Henry, Lionel, and Hadley Wickham. 2019. Purrr: Functional Programming Tools. https://CRAN.R-project.org/package=purrr.

Hothorn, Torsten. 2019. TH.data: TH’s Data Archive. https://CRAN.R-project.org/package=TH.data.

Hothorn, Torsten, Frank Bretz, and Peter Westfall. 2019. Multcomp: Simultaneous Inference in General Parametric Models. https://CRAN.R-project.org/package=multcomp.

Kuznetsova, Alexandra, Per Bruun Brockhoff, and Rune Haubo Bojesen Christensen. 2019. LmerTest: Tests in Linear Mixed Effects Models. https://CRAN.R-project.org/package=lmerTest.

Lenth, Russell. 2019. Emmeans: Estimated Marginal Means, Aka Least-Squares Means. https://CRAN.R-project.org/package=emmeans.

Müller, Kirill, and Hadley Wickham. 2019. Tibble: Simple Data Frames. https://CRAN.R-project.org/package=tibble.

R Core Team. 2020. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing. https://www.R-project.org/.

Ripley, Brian. 2019. MASS: Support Functions and Datasets for Venables and Ripley’s Mass. https://CRAN.R-project.org/package=MASS.

Therneau, Terry M. 2019. Survival: Survival Analysis. https://CRAN.R-project.org/package=survival.

Wickham, Hadley. 2017. Tidyverse: Easily Install and Load the ’Tidyverse’. https://CRAN.R-project.org/package=tidyverse.

———. 2019a. Forcats: Tools for Working with Categorical Variables (Factors). https://CRAN.R-project.org/package=forcats.

———. 2019b. Stringr: Simple, Consistent Wrappers for Common String Operations. https://CRAN.R-project.org/package=stringr.

Wickham, Hadley, and Lionel Henry. 2020. Tidyr: Tidy Messy Data. https://CRAN.R-project.org/package=tidyr.

Wickham, Hadley, Winston Chang, Lionel Henry, Thomas Lin Pedersen, Kohske Takahashi, Claus Wilke, Kara Woo, and Hiroaki Yutani. 2019. Ggplot2: Create Elegant Data Visualisations Using the Grammar of Graphics. https://CRAN.R-project.org/package=ggplot2.

Wickham, Hadley, Romain François, Lionel Henry, and Kirill Müller. 2019. Dplyr: A Grammar of Data Manipulation. https://CRAN.R-project.org/package=dplyr.

Wickham, Hadley, Jim Hester, and Romain Francois. 2018. Readr: Read Rectangular Text Data. https://CRAN.R-project.org/package=readr.